DRAFT¶
Ajitsuke Tamago - The art of boiling eggs¶
Steven Fowler (see more projects here)¶
TODO: Size this image correctly
Ajitsuke Tamago, aka Ramen eggs, are one of my favorite additions to a hearty ramen bowl. These eggs have a thick treacly yoke with a soft custard like white that has been marinated. The hardest part of making these eggs is the cooking to get the right egg consistencies.
The goal of this investigation is to look model the cooking step of making Ajitsuke Tamago eggs and get an understanding of most sensitive parameter that need to be controlled to get a repeatable product.
The chemistry of eggs¶
There are many publications on eggs, their chemistry and even modelling the boiling of eggs 1, 2. In summary, the composition and transformation of cooking an egg is as follows:
Based on Hen eggs, the main areas of an egg consists of 3 key parts:
- Shell - This is a porous structure mainly consisting of calcium carbonate and it accounts of 12% of an eggs mass.
- White (albumen) - Multi-layered opalescent substance mainly consisting of 2 proteins (12% ovotransferrin, 54% ovalbumin). It accounts for 67% of the egg's mass.
- Yoke - Contains fats, vitamins and < 50% of the egg's protein. It accounts for 33% of the egg's mass.
Other parts:
- Inner and outer membrane between shell and albumen.
- Chalazae - anchors the yoke in the centre of the egg.
- Vitelline membrane - between the yolk and albumen.
- Air cell - Generally forms at the larger end of the egg between the inner and outer membranes. It forms when the egg contracts after being laid.
How the parts of an egg cook (or not)¶
The phase change of an egg (yolk and white) is due to a coagulation of denatured proteins within the egg. The proteins denature at different temperature ranges, aggregate and form gel like structures. Note: Proteins can also be denatured by other mechanisms (i.e. pH changes, mechanical action) but these are out-of-scope of this investigation.
The phase changes of an egg are summarized in the below table (based on ref 3):
| Temperature (°C) | White | Yolk |
|---|---|---|
| 63 | Partially set | Runny thickening |
| 65 - 70 | Soft set gel | Runny treacle |
| 73 | Largely set | Soft gel |
| 77 | Tender solid | Hard boiled |
| 80 | Firm | Onset green discoloration |
| 90 | Over cooked | Dry, crumbly solid |
A general recipe for Ajitsuke Tamago is as follows4:
- Boil water (100°C)
- Place egg in water and cook for 7 minutes
- Immerse egg in ice water (0°C)
- Once cooled remove shells
- Marinade eggs for 24 hrs
Modelling requirements¶
Based on the above understanding, an Ajitsuke Tamago requires:
| ID_req | Requirement | Criteria |
|---|---|---|
| 1 | Yolk temperature shall not exceed 70°C to keep yolk runny | max(yolk temperature) =< 70°C |
| 2 | Yolk temperature shall exceed 65°C to thicken the yolk | mean(yolk temperature) > 65°C |
| 3 | White temperature shall not overly exceed 80°C to prevent onset of discoloration | mean(white temperature) =< 80°C |
Modelling Ajitsuke Tamago¶
The section provides an overview of the model used to simulate the Ajitsuke Tamago cooking. There are more detailed noted in the appendix. In summary, the following is covered:
- An overview of the maths deriving the model
- The creation of the mesh for the model
- Defining the egg shape and areas of interest
- Creating the mesh (diffusion coefficient array)
- The implementation of the model using the mesh
Overview of modelling math¶
To model the cooking of an egg, a general diffusion model was implemented to simulate the heat transfer in the system. A 2D diffusion type system in has the general form of: $$ \frac{\delta u}{\delta t} = \alpha (\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta i^2}) $$ where $u(x,i,t)$ is the function to be solved, $j$ and $i$ are cartesian coordinates and t is time. $\alpha$ is the diffustion coefficient and determines the rate change in $u$ in $t$.
Based on finite differences (see appendix for detailed notes), this can be expressed as: $$ u_i^{n+1} = u_{i,j}^n + F_x(u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n) + F_y(u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n) $$ where $$ F_x = \alpha \frac{\Delta t}{\Delta j^2}, F_y = \alpha \frac{\Delta t}{\Delta i^2} $$
[!Note] $F$ is the Fourier Number, a dimensionless value that represents the physical parameters of the diffusion problem in $\alpha$ and the discretized parameters $\Delta j$ and $\Delta t$.
Note the time domain is defined as: $$ t_n = n \Delta t, n = 0,...,N_t $$ with the space domain as: $$ x_i = i\Delta x, i = 0,...,N_x, $$ $$ y_j = j\Delta i, j = 0,...,N_y $$
Modelling Mesh¶
The standard mesh for a finite difference problem is a cartesian grid. Although there is the option to make more complex or dynamic meshes, if this is needed to capture the detail in an simulation then an FEA approach is probably more desirable. However, if we assume that the relatively high diffusion coefficient (≈ x10 greater than yolk and white) and a thickness of ≈0.3 mm of the shell has little its influence on the overall heat transfer, it can be ignored. So, the overall model consists of 2 parts: the white and the yolk.
Egg shape¶
The shape of the egg based on the model from Lian et al5, where the parametric form is: $$ y_p=asin\theta $$ $$ x_p=bcos\theta(1+c_1sin\theta+c_2sin^2\theta+c_3sin^3\theta) $$ where $y_p$ and $x_p$ are the coordinates of the egg profile (2D); $\theta$ is the eccentric angle; $a$ the half length; $b$ approximate half width; $c_1,c_2,c_3$ are parameters to be adjusted to fit an egg shape.
Show code
# Egg shape curve
# A parametric equation to generate an egg shape curve
# Reference: Comparison of egg-shape equations using relative curvature measures of nonlinearity, Lian et al., 2024
import numpy as np
import matplotlib.pyplot as plt
def egg_shape(half_width, half_height, t, c1, c2, c3 , x_axis_offset=0, y_axis_offset=0):
x = (half_width * np.cos(t)) * (1 + c1*np.sin(t) + c2*np.sin(t)**2 + c3*np.sin(t)**3) + x_axis_offset
i = half_height * np.sin(t) + y_axis_offset
return x, i
def find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , t, x_axis_offset, y_axis_offset):
egg_x, egg_y = egg_shape(half_width, half_height, t, c1, c2, c3, 0, 0)
err = np.atan2(point_y - y_axis_offset, point_x - x_axis_offset) - np.atan2(egg_y, egg_x)
# print("point: %f, %f - egg: %f, %f - err:, %f - theta: %f, %f"%(point_x, point_y, egg_x, egg_y, err, t, np.atan2(egg_y, egg_x)))
# print("theta guess: ", np.atan2(point_y - y_axis_offset, point_x - x_axis_offset))
return err, egg_x, egg_y
def find_theta_of_egg(point_x, point_y, half_width, half_height, c1 , c2 , c3 , x_axis_offset, y_axis_offset):
# Bisection method t = theta = atan2(i, x)
a = 0
b = np.pi
err_a = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , a, x_axis_offset, y_axis_offset)[0]
err_b = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , b, x_axis_offset, y_axis_offset)[0]
if np.sign(err_a) == np.sign(err_b):
# print("err_a: ", err_a, "a: ", a, "err_b: ", err_b, "b: ", b)
b = -np.pi
err_b = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , b, x_axis_offset, y_axis_offset)[0]
if np.sign(err_a) == np.sign(err_b):
# print("err_a: ", err_a, "a: ", a, "err_b: ", err_b, "b: ", b)
raise Exception("Bisection method requires a sign change over the interval.")
count = 0
while True:
m = (a + b) / 2
err_m, egg_x, egg_y = find_error_at_t(point_x, point_y, half_width, half_height, c1 , c2 , c3 , m, x_axis_offset, y_axis_offset)
if np.abs(err_m) < 1e-14:
# print("Converged in ", count, " iterations with error ", err_m)
return m, egg_x, egg_y
elif np.sign(err_a) == np.sign(err_m):
a = m
err_a = err_m
else:
b = m
err_b = err_m
count += 1
if count >= 1000:
print("Error: ", err_m)
raise Exception("Bisection method did not converge within 1000 iterations.")
def point_within_egg(point_x, point_y, half_width, half_height, c1 , c2 , c3 , x_axis_offset, y_axis_offset):
# Find a theta of the egg shape that matches the theta of the point
egg_theta, egg_x, egg_y = find_theta_of_egg(point_x=point_x, point_y=point_y, half_width=half_width, half_height=half_height,
c1=c1, c2=c2, c3=c3, x_axis_offset=x_axis_offset, y_axis_offset=y_axis_offset)
# Get egg shape X, Y that intersects egg_theta
# egg_x, egg_y = egg_shape(half_width, half_height, egg_theta, c1=c1, c2=c2, c3=c3, x_axis_offset=0, y_axis_offset=0)
# egg_x, egg_y = egg_shape(half_width, half_height, egg_theta, c1, c2, c3, 0, 0)
# print("theta point: ", np.atan2(point_y - y_axis_offset, point_x - x_axis_offset), ", theta of egg_cor: ", np.atan2(egg_y, egg_x), ", theta: ", egg_theta)
# print("egg x,i: ", egg_x, ", ", egg_y, " theta: ", np.atan2(egg_y, egg_x), " theta cor: ", egg_theta)
# Return if point is within egg shape
return ((point_x - x_axis_offset)**2 + (point_y - y_axis_offset)**2) <= (egg_x**2 + egg_y**2), egg_x, egg_y
def yoke_shape(x_axis_offset, y_axis_offset, a, b, t):
x = x_axis_offset + a*np.cos(t)
i = y_axis_offset + b*np.sin(t)
return x, i
###
#
# Testing egg shape functions and allocation of points
t = np.linspace(0, 2*np.pi, 100)
# Egg parameters
width = 15
height = 25
c1 = -0.2
c2 = 0.1
c3 = 0.0
x_shift = width
y_shift = height
# Yoke parameters
a = 10
b = 10
x_shift_yoke = x_shift
y_shift_yoke = y_shift - 10
x, y = egg_shape(width, height, t, c1=c1, c2=c2, c3=c3, x_axis_offset=x_shift, y_axis_offset=y_shift)
x_yoke, y_yoke = yoke_shape(x_shift_yoke, y_shift_yoke, a, b, t)
x_test = np.random.rand(1000) * width * 3 - 0.5 * width
y_test = np.random.rand(1000) * height * 3 - 0.5 * height
# x_test = np.ones(1) * 5
# y_test = np.ones(1) * 5
point_in_egg = np.zeros_like(x_test, dtype=bool)
point_in_york = np.zeros_like(x_test, dtype=bool)
x_egg_chk = np.zeros_like(x_test)
y_egg_chk = np.zeros_like(x_test)
for i in range(len(x_test)):
point_in_egg[i], x_egg_chk[i], y_egg_chk[i] = point_within_egg(x_test[i], y_test[i], width, height,
c1 , c2 , c3 ,
x_shift, y_shift)
# Need to fix this for an ellipse not circle
point_in_york[i] = ((x_test[i] - x_shift_yoke)**2 + (y_test[i] - y_shift_yoke)**2)**0.5 <= (a)
plt.plot(x_test[np.where(np.invert(point_in_egg))], y_test[np.where(np.invert(point_in_egg))], 'go') # Plot test points
plt.plot(x_test[np.where(point_in_egg != point_in_york)], y_test[np.where(point_in_egg != point_in_york)], 'ro') # Plot test points
plt.plot(x_test[np.where(point_in_york)], y_test[np.where(point_in_york)], 'yo') # Plot test points
plt.plot(x_egg_chk, y_egg_chk, 'bo') # Plot test points
plt.plot(x, y, 'b-')
plt.plot(x_yoke, y_yoke, 'b-')
plt.title('Egg Shape Curve - simulation using random points')
plt.xlabel('x')
plt.ylabel('i')
plt.axis('equal')
plt.style.use('seaborn-v0_8-whitegrid')
plt.show()
Figure 1. Illustration of egg shape fitting function.
The blue line shows the egg like shape that has been generated by the function. Randomly generated points are allocated to the specific area of the array/image: Yellow points are within the yolk; Red points are within the white and; green points are outside the egg shape. The blue points are those on the same vector as the randomly generated points (plotted from the origin).
The above plot shows the egg shape using fictional parameters to simulate an egg. The point in the plot are randomly generated and color-assigned depending on where are in the egg. This approach will be used to generate the mesh and thermal diffusion coefficient array for the modelling.
Define thermal diffusion coefficient array¶
The model will be defined by an array of varying diffusion coefficients that represent each area of the egg (white and yolk). Using the functions defined in the above section, the parameters of the egg fitting function were adjusted to fit measured data (see table below) and the yolk size was estimated based on the mass of the white and yolk (yolk being ≈33% of white + yolk area).
Three measures were sampled on medium sized supermarket eggs. The parameters of the egg shape equation were then adjusted to fit these measurements.
| Measurement | distance (mm) |
|---|---|
| Max width | 45.8 |
| Heigh of 'max width' | 25 |
| Max height | 55.8 |
The thermal properties of the egg were found in publications (see table below). However, as can be seen in the table the quoted magnitude of these properties varies considerably.
| Egg part | Reference | Density (Kg.cm-3) | Conductivity (W.m-1.K-1) | Capacity (J.Kg-1.K-1) |
|---|---|---|---|---|
| Yolk | 6 | 1028 | 0.54 | 3700 |
| Yolk | 7 | 1100 | 0.31 | 3700 |
| Yolk* | 8 | 1036 | 0.34 | 2093 |
| Yolk* | 8 | 1036 | 0.34 | 3266 |
| Yolk | 9 | N/A | 0.34 | 3180 |
| Yolk | 10 | 1032 | 0.34 | 2700 |
| White | 7 | 1134 | 0.48 | 3800 |
| White* | 8 | 1036 | 0.59 | 2093 |
| White* | 8 | 1036 | 0.59 | 3935 |
| White | 9 | N/A | 0.58 | 3180 |
| White | 10 | 1038 | 0.54 | 3700 |
For the purposes of this model, the values from reference [11] were chosen:
| Egg part | Density (Kg.cm-3) | Conductivity (W.m-1.K-1) | Capacity (J.Kg-1.K-1) | Diffusion (mm2.s-1) |
|---|---|---|---|---|
| Yolk | 1032 | 0.34 | 2700 | 0.12 |
| White | 1038 | 0.54 | 3700 | 0.14 |
Show code
import math
import numpy as np
def round_up_to_nearest_odd(value):
round_up_value = math.ceil(value)
if round_up_value % 2 == 0:
return int(round_up_value +1)
else:
return int(round_up_value)
# Make diffusion grid
# Model definition
egg_width = 45.8 # mm
egg_height = 55.8 # mm
egg_c1 = -0.1 # Shape parameter 1st order
egg_c2 = 0.1 # Shape parameter 2nd order
egg_c3 = 0.0 # Shape parameter 3rd order
white_alpha = 0.14060303077644118 # Egg white coefficient of diffusion
yoke_width = egg_width * 0.41**0.5 # mm
yoke_height = yoke_width # mm, assume round for now - can implement ellipse later if needed
yoke_center = -5 # mm, height from base off egg
yoke_alpha =0.1220212460522538 # Egg yoke coefficient of diffusion
# Mesh definition
dx = 0.5 # mm per node
dy = dx
# # mesh_width: odd value puts line of node at line of symmetry, +1 is for boundary conditions at 0 & n
# mesh_width = round_up_to_nearest_odd(egg_width / dx + 1)
# mesh_height = round_up_to_nearest_odd(egg_height / dy + 1)
mesh_X_mm = np.linspace(-egg_width/2 - 1, egg_width/2 + 1, int(np.ceil((egg_width + 2)/dx)))
mesh_Y_mm = np.linspace(-egg_height/2 - 1, egg_height/2 + 1, int(np.ceil((egg_height + 2)/dx)))
# mesh_X_mm = np.linspace(-egg_width/2 - 1, egg_width/2 + 1, (egg_width + 2)*dx)
# mesh_Y_mm = np.linspace(-egg_height/2 - 1, egg_height/2 + 1, (egg_width + 2)*dy)
# Diffusion grid
X, Y = np.meshgrid(mesh_X_mm, mesh_Y_mm, indexing='ij')
alpha_grid = np.zeros_like(X)
boundary_coord = []
white_coord = []
yolk_coord = []
for i in range(np.size(alpha_grid,0 )):
for j in range(np.size(alpha_grid, 1)):
# Check if point in yoke
point_in_yoke = (X[i, j]**2 + (Y[i, j] - yoke_center)**2) <= (yoke_width/2)**2
if point_in_yoke:
# alpha_grid[i, j] = yoke_alpha
yolk_coord.append((i, j))
else:
# Check if point in white
point_in_egg, e1, e2 = point_within_egg(X[i, j], Y[i, j], egg_width/2, egg_height/2,
egg_c1, egg_c2, egg_c3,
0, 0)
if point_in_egg:
# alpha_grid[i, j] = white_alpha
white_coord.append((i, j))
else:
boundary_coord.append((i, j))
boundary_coord = np.array(boundary_coord)
white_coord = np.array(white_coord)
yolk_coord = np.array(yolk_coord)
print(f"yolk to white ratio: {np.size(yolk_coord,0) / (np.size(white_coord,0) + np.size(yolk_coord,0)):.2f}")
i, j = egg_shape(egg_width/2, egg_height/2,
np.linspace(0, 2*np.pi, 100),
egg_c1, egg_c2, egg_c3,
0, 0)
alpha_grid[yolk_coord[:, 0], yolk_coord[:, 1]] = yoke_alpha
alpha_grid[white_coord[:, 0], white_coord[:, 1]] = white_alpha
alpha_grid[boundary_coord[:, 0], boundary_coord[:, 1]] = 0
# np.put(alpha_grid, np.ravel_multi_index(yolk_coord.T, alpha_grid.shape), diffusion_coefficient(25, yolk))
fig, ax = plt.subplots()
ax.set_aspect('equal')
scatter = ax.scatter(X, Y, c=alpha_grid, cmap='viridis', s=0.1)
fig.colorbar(scatter, ax=ax, label='Diffusion coefficient')
ax.plot(i, j)
ax.set_title('Image of thermal diffusion coefficient array')
ax.set_xlabel('X (mm)')
ax.set_ylabel('Y (mm)')
plt.show()
Figure 2. Illustration of diffusion coefficient array.
This plot shows the nodes of the cartesian mesh, where the colour of the node represents the magnitude of the diffusion coefficient. The blue line highlights the boundary of the egg. Note: point with 0 diffusion coefficient (dark purple) are the boundary nodes.
The above image illustrates the diffusion array for the model. Outside of the egg is the boundary nodes with a coefficient of 0, the yellow (≈0.11 mm2s-1) is the white and green (≈0.08 mm2s-1) is the yolk. The origin is at the centre of the egg-shape function and there is a 0.5 mm resolution.
Finite-difference model¶
A generic model for finite differences was drafted (see appendix) for solving general diffusion problem. However, for this investigation an diffusion coefficient array was needed to define the circular egg shape (without using polar coordinate system) and to look at the requirements the heating temperature needed to be dynamic. Hence, the below $ ajitsuke_egg_simulation $ function was written.
Criteria to investigate requirements
- Geometrically defined diffusion coefficients
- Dynamic heating (control of boundary conditions)
- Max temp monitoring (logic to drive heating temperature)
- Mean yolk monitoring (logic to drive heating temperature)
Model assumptions¶
| ID_ma | Assumption |
|---|---|
| 1 | Diffusion coefficients are constant throughout process (see appendix - this is not true) |
| 2 | The error in the finite-difference approach is within the noise of the system |
| 3 | The egg shape defined in the coefficient array is representative of a general egg (small samples numbers measured) |
| 4 | A uniform boundary condition will give representative results i.e. ignoring the dynamics of the boiling water |
| 5 | ...
Show code
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.cm as cm
import time as time
from IPython.display import HTML
import numpy as np
def ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid,
heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
boundary_T_cool, boundary_T_max, boundary_T_initial,
egg_T_initial, boundary_coord, yolk_coord, white_coord):
start_time = time.time()
nx = np.size(alpha_grid, 0)
ny = np.size(alpha_grid, 1)
# Check stability condition - using max alpha for worst case
if dt > (dx**2 * dy**2) / (2 * alpha_grid.max() * (dx**2 + dy**2)):
dt_max = (dx**2 * dy**2) / (2 * alpha_grid.max() * (dx**2 + dy**2))
alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)
# Set initial conditions
boundary_T_current = boundary_T_initial
u0 = np.empty_like(alpha_grid)
u0[boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
u0[white_coord[:, 0], white_coord[:, 1]] = egg_T_initial
u0[yolk_coord[:, 0], yolk_coord[:, 1]] = egg_T_initial
# Boundary conditions with time
cool = True
mean_yolk = None
mean_white = None
u = np.empty((nt, nx, ny)) # Initialize the solution array
u[0, :, :] = u0 # Set initial condition
for n in range(1, nt):
mean_yolk = np.mean(u[n-1, yolk_coord[:, 0], yolk_coord[:, 1]])
max_yolk = np.max(u[n-1, yolk_coord[:, 0], yolk_coord[:, 1]])
mean_white = np.mean(u[n-1, white_coord[:, 0], white_coord[:, 1]])
print(f"Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds\r", end='')
# Cut off heat source when requirement 2 is met (mean_yolk_temp_switch)
if ((mean_yolk >= mean_yolk_temp_switch) or (max_yolk >= max_yolk_temp_switch)) and (cool):
# u0[boundary_coord[:, 0], boundary_coord[:, 1]] = 0
# u[n, boundary_coord[:, 0], boundary_coord[:, 1]] = 0
boundary_T_current = boundary_T_cool
u[n-1, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
cool = False
print(f"Temperature changed! - Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds")
# Cap heating to boundary_T_max
if (boundary_T_current < boundary_T_max) & (cool):
boundary_T_current += heat_rate
u[n-1, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
# Boundary condition check
u[n, boundary_coord[:, 0], boundary_coord[:, 1]] = boundary_T_current
# Iterate over internal grid points
for i in range(1, nx-1):
for j in range(1, ny-1):
u[n, i, j] = (u[n-1, i, j] +
alpha_grid[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
(u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
end_time = time.time()
elapsed_time = end_time - start_time
print(f"Mean yolk: {mean_yolk:.2f}, Mean white: {mean_white:.2f} - after {dt*n:.2f} seconds")
print(f"Elapsed time for solution: {elapsed_time:.2f} seconds.\n Time simulated: {dt * nt:.2f} seconds.")
return u
Show code
dt = 0.035
nt = 35000 #8600
mean_yolk_temp_switch = 1000 # Too high to trigger
max_yolk_temp_switch = 1000 # Too high to trigger
boundary_T_cool = 0
boundary_T_max = 100
egg_T_initial = 20
###
heat_rate = 0 * dt
boundary_T_initial = 100
data_100 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid,
heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
boundary_T_cool, boundary_T_max, boundary_T_initial,
egg_T_initial, boundary_coord, yolk_coord, white_coord)
# Save for later access without needing to re-run simulation
np.save("egg_100_heating", data_100)
Show code
from matplotlib import collections as mc
###
# Option to load data rather than re-running simulation
###
# dt = 0.035
# data_100 = np.load("egg_100_heating.npy")
# TODO: Write plotting function based on below to tidy code and uniform plots
x = np.zeros(np.size(data_100, 0))
y_label = ["boundary T", "mean white T", "mean yolk T", "min yolk T", "max yolk T", "max_white_T"]
y_colour = ["g-", "y-", "b-", "b--", "b--", "w "]
y_100 = np.zeros((np.size(y_label), np.size(x)))
for n in range(np.size(x)):
x[n] = n * dt
y_100[0, n] = np.mean(data_100[n, boundary_coord[:, 0], boundary_coord[:, 1]])
y_100[1, n] = np.mean(data_100[n, white_coord[:, 0], white_coord[:, 1]])
y_100[2, n] = np.mean(data_100[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_100[3, n] = np.min(data_100[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_100[4, n] = np.max(data_100[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_100[5, n] = np.max(data_100[n, white_coord[:, 0], white_coord[:, 1]])
fig, ax = plt.subplots()
for n in range(np.size(y_label)):
ax.plot(x, y_100[n], y_colour[n], label=y_label[n])
ax.fill_between(x, y_100[3], y_100[4], facecolor='grey', alpha=0.5)
ax.hlines([63, 80], xmin=np.min(x), xmax=np.max(x), colors='k', linestyles='-.')
plt.text(0, 63, '63°C', verticalalignment='bottom', horizontalalignment='left', color='k')
plt.text(0, 80, '80°C', verticalalignment='bottom', horizontalalignment='left', color='k')
ax.set_title('Water at 100°C from t$_0$')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Temperature (°C)')
# plt.tight_layout()
fig.suptitle('Ahitsuke Tamago simulation:')
plt.show()
Figure 3. Simulation of boiling egg at 100°C for > 1200 seconds.
TODO: Update - Green is boundary condition (temperature of water cooking egg); Yellow is the mean temperature of the white; Blue is the mean (solid), max and min (dashed) yolk temperature. Grey area highlights temperature range of yolk. Black dash-dot lines are at: 63°C and 80°C.
Show code
# Useful numbers for discussion below...
time_reach_63 = x[np.argmax(y_100[4] >= 63)]
print("Data points of interest...\n")
print(f"\tx Yolk starts to heat (T > 21°C) at: {x[np.argmax(y_100[4] > 21)]:.2f} seconds")
print(f"\tx Yolk-white interface reaches 63°C: {time_reach_63:.2f} seconds ({time_reach_63/60:.1f} minutes)")
print(f"\tx Mean Yolk temperature when yolk-white interface at 63°C: {y_100[2, np.argmax(x >= time_reach_63)]:.2f}°C")
print(f"\tx Mean White temperature when yolk-white interface at 63°C: {y_100[1, np.argmax(x >= time_reach_63)]:.2f}°C")
print(f"\tx Max White temperature when yolk-white interface at 63°C: {y_100[5, np.argmax(x >= time_reach_63)]:.2f}°C")
From the above plot it can be seen that the yolk does not start to heat for 39 seconds, it them reaches 63°C at 7 minutes. This is a nice alignment as the recipe for cooking the Ajitsuke Tamago stated a 7 minutes cooking time and reference [6] [10] defines it's egg simulation complete when the yolk-white interface reaches 63°C. When the yolk-white interface reaches 63°C the mean white temperature is at 77°C, which is below at the 'tender' hardness (firm is at 80°C).
Heating rate to meet requirements¶
TODO: Check ref. 10 As discussed above, reference [6] [10] assumed that the egg is cooked when the yolk-white interface reached 63°C. Consequently, for the next simulations when this condition is met the boundary condition will be dropped to 0°C simulated the egg being submerged in ice water.
3 differing heating scenarios will be modelled:
- Scenario 1 - Starting with boiling water (reference)
- Scenario 2 - Starting at 20°C and heating to 100°C at a rate of 0.50°C.s-1
- Scenario 3 - Starting at 20°C and heating to 100°C at a rate of 0.25°C.s-1
This is to simulate the egg being present as the water is boiled.
Show code
dt = 0.035
nt = 35000 #8600
# dx
# dy
# alpha_grid # Defines dimensions of problem
heat_rate = 0.5 * dt
mean_yolk_temp_switch = 100 # Requirement 2
max_yolk_temp_switch = 63
boundary_T_cool = 0
boundary_T_max = 100
boundary_T_initial = 20
egg_T_initial = 20
# boundary_coord
# yolk_coord
# white_coord
###
data_050 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid,
heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
boundary_T_cool, boundary_T_max, boundary_T_initial,
egg_T_initial, boundary_coord, yolk_coord, white_coord)
np.save("egg_050_heating", data_050)
###
heat_rate = 0.25 * dt
data_025 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid,
heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
boundary_T_cool, boundary_T_max, boundary_T_initial,
egg_T_initial, boundary_coord, yolk_coord, white_coord)
np.save("egg_025_heating", data_025)
###
heat_rate = 0 * dt
boundary_T_initial = 100
data_000 = ajitsuke_egg_simulation(dt, nt, dx, dy, alpha_grid,
heat_rate, mean_yolk_temp_switch, max_yolk_temp_switch,
boundary_T_cool, boundary_T_max, boundary_T_initial,
egg_T_initial, boundary_coord, yolk_coord, white_coord)
np.save("egg_000_heating", data_000)
Show code
from matplotlib import collections as mc
###
# Option to load data rather than re-running simulation
###
# dt = 0.035
# data_050 = np.load("egg_050_heating.npy")
# data_025 = np.load("egg_025_heating.npy")
# data_000 = np.load("egg_000_heating.npy")
# TODO: Use plotting function
x = np.zeros(np.size(data_000, 0))
y_label = ["boundary T", "mean white T", "mean yolk T", "min yolk T", "max yolk T", "max_white_T"]
y_colour = ["g-", "y-", "b-", "b--", "b--", "w "]
y_000 = np.zeros((np.size(y_label), np.size(x)))
y_050 = np.zeros((np.size(y_label), np.size(x)))
y_025 = np.zeros((np.size(y_label), np.size(x)))
for n in range(np.size(x)):
x[n] = n * dt
y_000[0, n] = np.mean(data_000[n, boundary_coord[:, 0], boundary_coord[:, 1]])
y_000[1, n] = np.mean(data_000[n, white_coord[:, 0], white_coord[:, 1]])
y_000[2, n] = np.mean(data_000[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_000[3, n] = np.min(data_000[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_000[4, n] = np.max(data_000[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_000[5, n] = np.max(data_000[n, white_coord[:, 0], white_coord[:, 1]])
y_050[0, n] = np.mean(data_050[n, boundary_coord[:, 0], boundary_coord[:, 1]])
y_050[1, n] = np.mean(data_050[n, white_coord[:, 0], white_coord[:, 1]])
y_050[2, n] = np.mean(data_050[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_050[3, n] = np.min(data_050[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_050[4, n] = np.max(data_050[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_050[5, n] = np.max(data_050[n, white_coord[:, 0], white_coord[:, 1]])
y_025[0, n] = np.mean(data_025[n, boundary_coord[:, 0], boundary_coord[:, 1]])
y_025[1, n] = np.mean(data_025[n, white_coord[:, 0], white_coord[:, 1]])
y_025[2, n] = np.mean(data_025[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_025[3, n] = np.min(data_025[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_025[4, n] = np.max(data_025[n, yolk_coord[:, 0], yolk_coord[:, 1]])
y_025[5, n] = np.max(data_025[n, white_coord[:, 0], white_coord[:, 1]])
fig, axs = plt.subplots(1, 3, figsize=(12, 5), sharey=True)
for n in range(np.size(y_label)):
axs[0].plot(x, y_000[n], y_colour[n], label=y_label[n])
axs[1].plot(x, y_050[n], y_colour[n], label=y_label[n])
axs[2].plot(x, y_025[n], y_colour[n], label=y_label[n])
axs[0].fill_between(x, y_000[3], y_000[4], facecolor='grey', alpha=0.5)
axs[1].fill_between(x, y_050[3], y_050[4], facecolor='grey', alpha=0.5)
axs[2].fill_between(x, y_025[3], y_025[4], facecolor='grey', alpha=0.5)
axs[0].hlines([63, 80], xmin=np.min(x), xmax=np.max(x), colors='k', linestyles='-.')
axs[1].hlines([63, 80], xmin=np.min(x), xmax=np.max(x), colors='k', linestyles='-.')
axs[2].hlines([63, 80], xmin=np.min(x), xmax=np.max(x), colors='k', linestyles='-.')
axs[0].text(0, 63, '63°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[0].text(0, 80, '80°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[1].text(0, 63, '63°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[1].text(0, 80, '80°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[2].text(0, 63, '63°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[2].text(0, 80, '80°C', verticalalignment='bottom', horizontalalignment='left', color='k')
axs[0].set_title('Water at 100°C from t$_0$')
axs[0].set_xlabel('Time (s)')
axs[0].set_ylabel('Temperature (°C)')
axs[1].set_title('Water heated at 0.5°C/s, t$_0$ 25°C')
axs[1].set_xlabel('Time (s)')
axs[1].set_ylabel('Temperature (°C)')
axs[2].set_title('Water heated at 0.25°C/s, t$_0$ 25°C')
axs[2].set_xlabel('Time (s)')
axs[2].set_ylabel('Temperature (°C)')
# plt.tight_layout()
fig.suptitle('Ahitsuke Tamago simulation:')
plt.show()
Figure 4. Ahitsuke Tamago cooking simulation with differing boundary conditions.
Update - Green is boundary condition (temperature of water cooking egg); Yellow is the mean temperature of the white; Blue is the mean (solid), max and min (dashed) yolk temperature. Grey area highlights temperature range of yolk. Black dash-dot lines are at 65°C and 70°C, the limits for the yolk.
Show code
# Useful numbers for discussion below...
print("Data points of interest...\n")
print(f"Maximum yolk temperature for: \n\tScenario 1 -> {np.max(y_000[4]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[4]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[4]):.2f}°C.\n")
print(f"Maximum white temperature for: \n\tScenario 1 -> {np.max(y_000[5]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[5]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[5]):.2f}°C.\n")
print(f"Maximum mean white temperature reached for: \n\tScenario 1 -> {np.max(y_000[1]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[1]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[1]):.2f}°C.\n")
print(f"Minimum egg (yolk) temperature reached for: \n\tScenario 1 -> {np.max(y_000[3]):.2f}°C; \n\tScenario 2 -> {np.max(y_050[3]):.2f}°C; \n\tScenario 3 -> {np.max(y_025[3]):.2f}°C.\n")
print(f"Time of 63°C inflection point: \n\tScenario 1 -> {x[np.argmax(y_000[0] == 0)]/60:.1f} min; \n\tScenario 2 -> {x[np.argmax(y_050[0] == 0)]/60:.1f} min; \n\tScenario 3 -> {x[np.argmax(y_025[0] == 0)]/60:.1f} min.\n")
Looking at the plot and reinforced by the 'useful numbers' above, apart from the initial heating of the egg the response and temperature at and after the 63°C inflection are all the same for the 3 scenarios. Consequently, the only difference the length of time until the 63°C inflection points which are at 7.0, 8.4 & 9.9 minutes for Scenario 1, 2, & 3 respectively.
Requirements¶
Considering the above plots in context of the requirements:
Yolk temperature shall not exceed 70°C to keep yolk runny
- This was met for all scenarios as the maximum yolk temperature was ≈ 64°C.
Yolk temperature shall exceed 65°C to thicken the yolk
- This was not met for all scenarios as the maximum temperature was ≈ 64°C.
White temperature shall not overly exceed 80°C to prevent onset of discoloration
- This was met for all scenarios as the maximum mean white temperature was ≈ 77°C.
In summary, all scenarios produce the same egg. However, it could be argued that scenario 1 is less prone to error as you don't need to factor the rate of water heating into the cooking.
Heat profile¶
TODO: Complete Staring in boiling water works and should deliver the required results (like presented in the recipe) Add summary of egg images and video
Show code
# Profile of egg at 63°C inflection point
index_reach_63 = np.argmax(y_000[4] >= 63)
crosssection_point = 50
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.text(5, crosssection_point, 'Cross-section', verticalalignment='bottom', horizontalalignment='left', color='b')
ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')
im1 = ax1.imshow(data_000[index_reach_63,].T, cmap='hot', vmin=0, vmax=100, origin='lower') # Transpose needed as array is [i,j] not [x,y]
# ax1.legend()
ax1.set_xlabel('Size (nodes)')
ax1.set_ylabel('Size (nodes)')
ax1.set_title('Heat map of egg (°C)')
fig.colorbar(im1, ax=ax1)
x_data = np.zeros(np.size(data_000,1))
y_data = np.zeros_like(x_data)
for i in range(np.size(x_data)):
x_data[i] = i * dx
y_data[i] = data_000[index_reach_63, i, crosssection_point]
ax2.plot(x_data, y_data, 'b-')
ax2.set_title('Cross section')
ax2.set_xlabel('Size (mm)')
ax2.set_ylabel('Temperaure (°C)')
# ax2.axvline((np.max(x_data) - yoke_width)/2, linestyle='-.')
# ax2.axvline((np.max(x_data) + yoke_width)/2, linestyle='-.')
# ax2.text((np.max(x_data) - yoke_width)/2 + 0.5, 95, 'Yoke-white\ninterface', verticalalignment='center', horizontalalignment='left', color='k')
fig.suptitle('Profile of egg at 63°C inflection point')
plt.tight_layout()
plt.show()
Figure 5. Heat map and cross-section of egg at 63°C inflection point.
Heat map shows 0°C boundary condition simulating egg's submersion into ice water. The profile shows the center of the yolk ≈ 30°C.
Show code
ani = None
count = 0
images = []
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')
ax2.set_ylim(-5, 105)
frames = np.arange(0, np.size(data_000, 0), 350)
x_data = frames * dt
y_data = np.empty((np.size(y_label), np.size(x_data)))
y_data[:, :] = np.nan
for n in range(np.size(y_data, 1)):
img = ax1.imshow(data_000[frames[n]].T, cmap='hot', vmin=0, vmax=100, animated=True, origin='lower')
line_art = None
# x_data[n] = frames[n] * dt
# for i in range(np.size(y_data, 0)):
for i in range(np.size(y_label)):
y_data[i, n] = y_000[i, frames[n]]
line_art0, = ax2.plot(x_data, y_data[0], y_colour[0], label=y_label[0], animated=True)
line_art1, = ax2.plot(x_data, y_data[1], y_colour[1], label=y_label[1], animated=True)
line_art2, = ax2.plot(x_data, y_data[2], y_colour[2], label=y_label[2], animated=True)
line_art3, = ax2.plot(x_data, y_data[3], y_colour[3], label=y_label[3], animated=True)
line_art4, = ax2.plot(x_data, y_data[4], y_colour[4], label=y_label[4], animated=True)
fill_art = ax2.fill_between(x_data, y_data[3], y_data[4], facecolor='grey', alpha=0.5, animated=True)
images.append([img, line_art0, line_art1, line_art2, line_art3, line_art4, fill_art])
count += 1
print(f"count: {count}\r", end='')
# plt.tight_layout()
ani = animation.ArtistAnimation(fig, images, interval=50, blit=True, repeat_delay=2000)
plt.close()
HTML(ani.to_jshtml())
Show code
ani = None
count = 0
images = []
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.axhline(crosssection_point, linestyle='--', label='plot line', color='blue')
for n in range(0, np.size(data_000,0), 350):
img = ax1.imshow(data_000[n].T, cmap='hot', vmin=0, vmax=100, animated=True, origin='lower')
x_data = np.zeros(np.size(data_000,1))
y_data = np.zeros_like(x_data)
for i in range(np.size(x_data)):
x_data[i] = i * dx
y_data[i] = data_000[n, i, crosssection_point]
line_art, = ax2.plot(x_data, y_data, 'b-', animated=True)
# ax2.axvline(20, linestyle='--', animated=True)
# ax2.axvline(40, linestyle='--', animated=True)
images.append([img, line_art])
count += 1
print(f"count: {count}\r", end='')
# plt.tight_layout()
ani = animation.ArtistAnimation(fig, images, interval=50, blit=True, repeat_delay=1000)
plt.close()
HTML(ani.to_jshtml())
References¶
[1]: https://www.sciencedirect.com/science/article/pii/S2665927124001989 - Modeling of cooking and phase change of egg white using computational fluid dynamics. R E. Sánchez-García et al. 2024
[2]: https://www.nature.com/articles/s44172-024-00334-w - Periodic cooking of eggs. E D. Lorenzo et al. 2025
[3]: https://math.arizona.edu/~gabitov/teaching/141/math_485/Final_Presentations/Cooking_Final_Presentation.pdf - Numerical Modeling of Diffusion and Phase Transitions in Heterogeneous Media. H. Jeffrey.
[4]: https://tokusen.store/en/blogs/recipes/ajitsuke-tamago-oeufs-marines-a-la-japonaise?srsltid=AfmBOop-R7FBSF8x4BBuf1P6eeOMCkhl4dBXGh3gvz_CM9kgHepp4Yd1 - Japanese Ramen Eggs (Ajitsuke Tamago). Tokusen [website]. 2025.
[4]: https://khymos.org/2009/04/09/towards-the-perfect-soft-boiled-egg/ - Towards the perfect soft boiled egg, M. Lersch. 2009
[5]: https://www.sciencedirect.com/science/article/pii/S0032579124006485 - Comparison of egg-shape equations using relative curvature measures of nonlinearity. M. Lain et al. 2024
[6]: https://newton.ex.ac.uk/teaching/CDHW/egg/ -
[7]: https://www.sciencedirect.com/science/article/pii/S2665927124001989 -
[8]: https://www.researchgate.net/figure/Thermal-properties-of-egg-components_tbl1_256804375 -
[9]: https://www.engineeringtoolbox.com/food-thermal-conductivity-d_2177.html -
[10]: https://newton.ex.ac.uk/teaching/CDHW/egg/CW061201-1.pdf -
Appendix¶
Notes from Lecture videos¶
Definition of a derivative $$ \left. \frac{df(x)}{dx} \right|_{x=0} = f'(a)=\lim_{x \to 0}\frac{f(x)-f(a)}{x-a} $$
Forward difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i} $$
Backward difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_i)-f(x_{i-1})}{x_i-x_{i-1}} $$
Central difference: $$ \left. \frac{df(x)}{dx} \right|_{x_i} \approx \frac{f(x_{i+1})-f(x_{i-1})}{x_{i+1}-x_{i-1}} $$
Taylor series¶
- Assume all nodes/points are equally spaced i.e. h is constant
$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(x_i) + \frac{1}{3!}h^3f'''(x_i) + \frac{1}{4!}h^4f''''(x_i) + ... $$
$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(x_i) + ... + \left. \frac{1}{n!}h^n \frac{d^nf}{x^n} \right|_{x_i} + R_n \xi $$
where $\xi$ is an unknown location between $x_i \ge \xi \ge x_{i+1}$ and we cannot calculate $\xi$.
Using the first 2 terms of the Taylor's Series for $f(x_{i+1})$
$$ f(x_{i+1}) = f(x_i) + hf'(x_i) + \frac{1}{2!}h^2f''(\xi) $$
Solve for $f'(x_{i+1})$: $$ f'(x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{h} - \frac{1}{2!}hf''(\xi) = \frac{f(x_{i+1})-f(x_i)}{h} - error $$
Note: As the value of $h$ gets smaller the $error$ is smaller: $$ error = \frac{f''(\xi)}{2}h $$ This is also a first order method (order of h given to the truncation error)
General 1D diffusion problem¶
Diffusion type systems have a general form of: $$ \frac{\delta u}{\delta t} = \alpha \frac{\delta^2 u}{\delta x^2} + f $$
The forward Euler model [^1], the space domain is defined with the following nodes: $$ x_i = i\Delta x, i = 0,...,N_x $$ and the time domain is defines as: $$ t_n = n \Delta t, n = 0,...,N_t $$
So, considering the general diffusion system notes above, $u(x_i, t_n)$ can be denoted as $u^n_i$. Thus, the general diffusion system becomes: $$ \frac{\delta}{\delta t}u(x_i, t_n) = \alpha \frac{\delta^2}{\delta x^2} u(x_i, t_n) $$
The solution to a 2^nd order finite difference approximation can be derived by combining 2 Taylor series, this results in a forward difference in time and a central difference in space: $$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n +u_{i-1}^n}{\delta x^2} + Err_{trunc}(dx^2) $$
Where the truncation error ($Err_{trunc}$) can be seen to be 2^nd order. Rearranging this with respect to the unknown $u_i^{n+1}$: $$ u_i^{n+1} = u_i^n + F(u_{i+1}^n - 2u_i^n + u_{i-1}^n) +\Delta t.Err_{trunc}(dx^2) $$ where: $$ F = \alpha \frac{\Delta t}{\Delta x^2} $$
[!Note] F is the Fourier Number, a dimensionless value that represents the physical parameters of the diffusion problem in $\alpha$ and the discretised parameters $\Delta x$ and $\Delta t$.
[^1]: Finite difference methods for diffusion processes. H. P. Langtangen, S. Linge. 2016 July 14
General 2D diffusion problem¶
Diffusion type systems in 2D have a general form of: $$ \frac{\delta u}{\delta t} = \alpha (\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta i^2}) $$ For the 2D problem we add another dimension: $$ y_j = j\Delta i, j = 0,...,N_y $$ So our 2D diffusion system, $u(x_i, y_j, t_n)$ can be denoted as $u^n_{i,j}$, thus: $$ u_i^{n+1} = u_{i,j}^n + F_x(u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n) + F_y(u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n) $$ where: $$ F_x = \alpha \frac{\Delta t}{\Delta x^2}, F_y = \alpha \frac{\Delta t}{\Delta i^2} $$ If we assume $\Delta x = \Delta i$, the 2D solution can be simplified to: $$ u_i^{n+1} = u_{i,j}^n + F(u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n) $$ where: $$ F = \alpha \frac{\Delta t}{\Delta x^2} $$
Simulation¶
...
Show code
""""
Function for solving diffusion system in 2D using finite difference method.
This function has the addition of:
- Variable alpha values for diffusion coefficient
"""
import numpy as np
def solve_diffusion_2d_vAlpha(nx, ny, dx, dy, dt, nt, alpha, u0):
"""
Solve the 2D diffusion equation using finite difference method.
Parameters:
nx : int
Number of grid points in x direction.
ny : int
Number of grid points in i direction.
dx : float
Grid spacing in x direction.
dy : float
Grid spacing in i direction.
dt : float
Time step size.
nt : int
Number of time steps to simulate.
alpha : 2D array
Diffusion coefficient.
u0 : 2D array
Initial condition for the diffusion equation.
Returns:
u : 3D array (t, x, i)
All solutions from t = 0 to nt time steps.
"""
# Check stability condition - using max alpha for worst case
if dt > (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2)):
dt_max = (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2))
alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)
u = np.zeros((nt, nx, ny)) # Initialize the solution array
u[0, :, :] = u0 # Set initial condition
for n in range(1, nt):
# Reset boundary conditions at each time step
u[n, 0, :] = u0[0, :] # Left boundary
u[n, -1, :] = u0[-1, :] # Right boundary
u[n, :, 0] = u0[:, 0] # Bottom boundary
u[n, :, -1] = u0[:, -1] # Top boundary
# Iterate over internal grid points
for i in range(1, nx-1):
for j in range(1, ny-1):
u[n, i, j] = (u[n-1, i, j] +
alpha[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
(u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
return u
Variable diffusion coefficient¶
As mentioned in the above, the properties of the egg (white & yolk) change with temperature, this including their effective thermal diffusion coefficient.
Show code
# Check diffusion coefficient
# https://newton.ex.ac.uk/teaching/CDHW/egg/
density_1 = 1.028 # g / cm^3
conduc_1 = 5.4e-3 # W /cm / K
capaci_1 = 3.7 # J /g /K
diffCo_1 = conduc_1 / (density_1 * capaci_1) # cm^2 / s
diffCo_1 = diffCo_1 * 100 # mm^2 /s
print(diffCo_1)
# https://www.sciencedirect.com/science/article/pii/S2665927124001989
print(1e6 * 0.31 / (1100 * 3700)) # yolk
print(1e6 * 0.48 / (1134 * 3800)) # white
print(1e6 * 0.6 / (1023 * 2600))
print(1e6 * 0.4 / (1143 * 3700))
# https://www.researchgate.net/figure/Thermal-properties-of-egg-components_tbl1_256804375
print(1e6 * 0.3370 / (1036 * 2093)) # yolk
print(1e6 * 0.3370 / (1036 * 3266)) # yolk
print(1e6 * 0.59 / (1036 * 2093)) # white
print(1e6 * 0.59 / (1036 * 3935)) # white
# https://www.engineeringtoolbox.com/food-thermal-conductivity-d_2177.html
# conduct 0.58 white, 0.34 yolk
# cap 3180
print(1e6 * 0.34 / (1036 * 3266))
print(1e6 * 0.58 / (1134 * 3000))
# https://newton.ex.ac.uk/teaching/CDHW/egg/CW061201-1.pdf
print(1e6 * 0.34 / (1032 * 2700)) # Yolk
print(1e6 * 0.54 / (1038 * 3700)) # White
Show code
# Modelling thermal coefficient as function of temperature (coagulation)
import numpy as np
import matplotlib.pyplot as plt
def thermal_diffusion(conductivity, density, capacity):
"""Return thermal diffusion coefficient
Args:
conductivity (float): thermal conductivity (k, W m-1 K-1)
density (float): density ('rho', kg m-3)
capacity (float): heat capacity (c, J kg-1 K-1)
Returns:
float: thermal diffusion ('alpha', m2 / s)
"""
return conductivity / (density * capacity)
def logistic(L, k, x, x0, y0):
"""Return logistic data
Args:
L (float): supremum (end value)
k (float): growth rate (logistic)
x (array, float): x value(s)
x0 (float): mid-point
y0 (float): i intercept
"""
return (L - y0) / (1 + np.exp(-k*(x - x0))) + y0
def diffusion_coefficient_mm2(T, material, width=0.001):
L = thermal_diffusion(material["conductivity"][0], material["density"][0], material["capacity"][0])
k = np.log((1-width)/width) / material["interval_K"]
x0 = material["chang_T"] + material["interval_K"]/2
y0 = thermal_diffusion(material["conductivity"][1], material["density"][1], material["capacity"][1])
return (1e6 * logistic(L, k, T, x0, y0))
# Data from https://www.sciencedirect.com/science/article/pii/S2665927124001989
yolk = {
"density": [1100, 1032], # kg m-3, [liquid, solid]
"conductivity": [0.31, 0.34], # W m-1 K-1
"capacity": [3700, 2700], # J kg-1 K-1
"chang_T": 343.15, # K
"interval_K": 3, # K
"latent":12.2 # KJ kg-1
}
white = {
"density": [1134, 1100], # kg m-3, [liquid, solid]
"conductivity": [0.48, 0.40], # W m-1 K-1
"capacity": [3800, 3700], # J kg-1 K-1
"chang_T": 336.15, # K
"interval_K": 3, # K
"latent":12.2 # KJ kg-1
}
shell = {
"density": [2300, 2300], # kg m-3, [solid, solid] <- this is N/A
"conductivity": [2.25, 2.25], # W m-1 K-1
"capacity": [888, 888], # J kg-1 K-1
"chang_T": 1, # K <- N/A
"interval_K": 1, # K <- N/A
"latent": 0 # KJ kg-1 <- N/A
}
# yolk["diffusion"] = [thermal_diffusion(yolk["conductivity"][0], yolk["density"][0], yolk["capacity"][0])]
# yolk["diffusion"].append(thermal_diffusion(yolk["conductivity"][1], yolk["density"][1], yolk["capacity"][1]))
# white["diffusion"] = [thermal_diffusion(white["conductivity"][0], white["density"][0], white["capacity"][0])]
# white["diffusion"].append(thermal_diffusion(white["conductivity"][1], white["density"][1], white["capacity"][1]))
# shell["diffusion"] = [thermal_diffusion(shell["conductivity"][0], shell["density"][0], shell["capacity"][0])]
# shell["diffusion"].append(thermal_diffusion(shell["conductivity"][1], shell["density"][1], shell["capacity"][1]))
# Plot diffusion coefficeint vs temperature, assume sigmodal change in coefficient
T = np.linspace(273.15, 273.15 + 100, 101)
# Assume a logistic change (fairly linear), with curve width (3K) fitting the 0.1% to 99.9% window
# yolk_diffusion = logistic(yolk["diffusion"][1], np.log(0.999/0.001) / yolk["interval_K"], T, yolk["chang_T"] + yolk["interval_K"]/2, yolk["diffusion"][0])
# white_diffusion = logistic(white["diffusion"][1], np.log(0.999/0.001) / white["interval_K"], T, white["chang_T"] + white["interval_K"]/2, white["diffusion"][0])
# shell_diffusion = logistic(shell["diffusion"][1], np.log(0.999/0.001) / shell["interval_K"], T, shell["chang_T"] + shell["interval_K"]/2, shell["diffusion"][0])
fig, ax = plt.subplots()
ax.plot(T, diffusion_coefficient_mm2(T, yolk), label="yolk", color='red')
ax.plot(T, diffusion_coefficient_mm2(T, white), label="white", color='blue')
plt.plot(T, diffusion_coefficient_mm2(T, shell), label="shell")
ax.legend()
ax.set_title("Change in diffusion coefficient with temperature")
ax.set_xlabel("Temperature (K)")
ax.set_ylabel("Thermal diffusion coefficient (mm$^2$.s$^{-1}$)")
plt.style.use('seaborn-v0_8-whitegrid')
plt.show()
The above plot shows how the diffusion coefficient changes as temperature increases. It can be noted that:
- The shell had a constant thermal diffusion coefficient
- The white's thermal diffusion coefficient decreases at the 63℃ transition temperature
- Note: This is assumed to be the transition due to the denaturing of the ovotransferrin protein as the later denaturing is not until 80℃.
- Therefore further change in diffusion coefficient > 80℃ are not covered in the data for this model.
- The yolk's thermal diffusion coefficient increases at the 70℃ transition temperature.
Show code
""""
Function for solving diffusion system in 2D using finite difference method.
This function has the addition of:
- Variable alpha values for diffusion coefficient
"""
import numpy as np
def solve_diffusion_2d_vAlpha(nx, ny, dx, dy, dt, nt, alpha, u0):
"""
Solve the 2D diffusion equation using finite difference method.
Parameters:
nx : int
Number of grid points in x direction.
ny : int
Number of grid points in i direction.
dx : float
Grid spacing in x direction.
dy : float
Grid spacing in i direction.
dt : float
Time step size.
nt : int
Number of time steps to simulate.
alpha : 2D array
Diffusion coefficient.
u0 : 2D array
Initial condition for the diffusion equation.
Returns:
u : 3D array (t, x, i)
All solutions from t = 0 to nt time steps.
"""
# Check stability condition - using max alpha for worst case
if dt > (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2)):
dt_max = (dx**2 * dy**2) / (2 * alpha.max() * (dx**2 + dy**2))
alpha_max = (dx**2 * dy**2) / (2 * dt * (dx**2 + dy**2))
raise ValueError("Time step size dt is too large for stability. Maximum dt: ", dt_max, " or reduce alpha to: ", alpha_max)
u = np.zeros((nt, nx, ny)) # Initialize the solution array
u[0, :, :] = u0 # Set initial condition
for n in range(1, nt):
# Reset boundary conditions at each time step
u[n, 0, :] = u0[0, :] # Left boundary
u[n, -1, :] = u0[-1, :] # Right boundary
u[n, :, 0] = u0[:, 0] # Bottom boundary
u[n, :, -1] = u0[:, -1] # Top boundary
# Iterate over internal grid points
for i in range(1, nx-1):
for j in range(1, ny-1):
u[n, i, j] = (u[n-1, i, j] +
alpha[i, j] * dt * ((u[n-1, i+1, j] - 2*u[n-1, i, j] + u[n-1, i-1, j]) / dx**2 +
(u[n-1, i, j+1] - 2*u[n-1, i, j] + u[n-1, i, j-1]) / dy**2))
return u
Show code
# Test alpha_grid
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.cm as cm
import time as time
from IPython.display import HTML
import numpy as np
initial_condition = np.ones_like(alpha_grid) * 25 # All start at 25 degrees
initial_condition[np.where(alpha_grid == 0)] = 100
# dt = 0.00625
nt = 480
start_time = time.time()
test = solve_diffusion_2d_vAlpha(np.size(alpha_grid, 1), np.size(alpha_grid, 0),
dx, dy, dt, nt, alpha_grid.T, initial_condition.T)
end_time = time.time()
elapsed_time = end_time - start_time
print(f"Elapsed time for solution: {elapsed_time:.2f} seconds.\n Time simulated: {dt * nt} seconds.")
# fig, (ax1, ax2) = plt.subplots(1, 2)
# images = []
# ax1.axhline(45, linestyle=':', label='plot line')
# im1 = ax1.imshow(test[-1].T, cmap='hot', vmin=25, vmax=100, origin='lower')
# ax1.legend()
# fig.colorbar(im1, ax=ax1)
# y_data = np.zeros(np.size(test,1))
# x_data = np.zeros_like(y_data)
# for j in range(np.size(test,1)):
# x_data[j] = j * dx
# y_data[j] = test[-1,j,45]
# ax2.plot(x_data, y_data, 'b-', animated=True)
# ax2.axvline(20, linestyle='--', animated=True)
# ax2.axvline(40, linestyle='--', animated=True)
# plt.show()
# plt.tight_layout()
# for i in range(100):
# ax1.axhline(45, linestyle=':', label='test')
# img = ax1.imshow(test[i*3].T, cmap='tab20', vmin=25, vmax=100, animated=True, origin='lower')
# y_data = np.zeros(np.size(test,1))
# x_data = np.zeros_like(y_data)
# for j in range(np.size(test,1)):
# x_data[j] = j
# y_data[j] = test[i,j,45]
# line_art, = ax2.plot(x_data, y_data, 'b-', animated=True)
# ax2.axvline(20, linestyle='--', animated=True)
# ax2.axvline(40, linestyle='--', animated=True)
# images.append([img, line_art])
# plt.tight_layout()
# ani = animation.ArtistAnimation(fig, images, interval=50, blit=True, repeat_delay=1000)
# HTML(ani.to_jshtml())
Stability¶
Add notes $$ \Delta t \leq \frac{\Delta x^2}{2\alpha} $$